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Abstract: 


The  incompressible  laminar  boundary  layer 
equations  are  solved  for  flow  about  a  sphere 
and  an  ellipsoid  using  a  finite  difference 
scheme.  The  effects  of  vortex  stretching  on 
the  momentum  and  vorticity  equations  are 
studied. 
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CHAPTER  t 


INTRODUCTION 

The  primary  purpose  of  this  study  is  to  develop  a  rapid  arid 
accurate  method  for  solving  the  incompressible  laminar  boundary 
layer  equations  for  axi symmetri c  flow  using  a  digital  computer. 

In  order  to  obtain  some  qualitative  insight  into  the  trends  exhibited 
by  model  calculations,  the  effects  of  vortex  stretching  on  the 
momentum  and  vorticity  equations  are  studied. 

There  are  many  methods  available  for  obtaining  solutions  to 
the  boundary  layer  equations;  a  fact  easily  verified  by  consulting 
a  standard  reference  such  as  Meksyn  (I96I),  Rosenhead  (1963)  or 
Schlichting  (I960).  The  majority  of  these  techniques  can  be 
classified  as  follows:  (1)  momentum  integral  methods,  (2)  correla¬ 
tion  methods,  (3)  infinite  series  method  and  (m)  finite  difference 
procedures.  In  general,  the  first  two  techniques  are  approximate 
procedures  in  the  sense  that  the  original  equations  have  been 
compromised.  The  infinite  series  and  finite  difference  methods, 
on  the  other  hand,  are  considered  to  be  theoretically  exact  in  the 
limit,  that  is,  for  an  'nfinite  number  of  terms  and  vanishingly 
small  step  lengths,  respectively.  Since  it  is  impossible  to 
discuss  all  of  the  individual  techniques  in  each  category,  only  a 
few  representative  examples  are  discussed  be  1 ow . 


The  Pohlhausen  and  Thwaites  methods  are  well-known 
examples  of  the  integral  and  correlation  techniques.  In  the  first 
of  these,  a  polynomial  expansion  for  the  velocity  is  assumed  and 
then  the  momentum  integral  aquation  is  solved.  Thwaites1  procedure 
is  based  on  an  empirical  correlation  of  existing  exact  and  approxi¬ 
mate  flow  solutions.  Using  this  known  data,  Thwaites  tabulated 
various  boundary  layer  shape  parameters  from  which  approximate 
values  for  the  momentum  and  displacement  thicknesses  can  be  readily 
determi ned. 
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The  method  described  by  Frossling  is  typical  of  the 
infinite  series  type  of  solution.  Here,  the  body  contour  and  the 
potential  flow  are  expressed  as  power  series  in  x  (distance 
measured  along  the  body).  The  stream  function  is  also  expanded 
in  an  infinite  series  in  xwith  coefficients  depending  on  the  waP 
distance  y.  These  series  are  then  substituted  into  the  boundary 
layer  equations,  written  in  terms  of  the  stream  function  and, 
Consequently,  an  order  set  of  ordinary  differential  equations 
results  for  the  unknown  coefficients. 

Finally,  the  Hartree-Womers  ley"*  scheme  is  considered  as  an 
illustration  of  the  finite  difference  type.  Briefly,  this  technique 
utilizes  a  simple  transformation  of  the  momentum  equation  in  which 
x,  y  remain  essentially  as  independent  variables.  The  region  of 
the  boundary  layer  is  divided  into  vertical  strips  and  the  trans¬ 
formed  momertum  equation  is  solved  numerically  at  each  step. 


3 


The  method  chosen  for  this  paper  is  based  on  a  finite 
difference  procedure  developed  by  A.  M.  0.  Smith  and  Darwin  U. 
Clutter^  which  is  actually  an  extension  of  the  Hart ree-Womers lay 
scheme.  The  most  important  reasons  for  its  selection  are 
(I)  it  is  more  accurate  for  a  wider  variety  of  problems  than  the 
approximate  integral  or  correlation  techniques,  (2)  there  is 
considerably  less  computational  work  involved  compared  to  the 
infinite  series  approach  when  high  accuracy  is  required,  (3)  it 
exhibits  a  greater  degree  of  numerical  stability  than  other  finite 
difference  methods,  and  (4)  the  equations  on  which  this  procedure 
is  based  are  free  of  singularities.  In  the  adaptation  of  Smith's 
and  Clutter's  method,  the  "starting"  procedure  (see  Reference  13 
and  Page  20  of  this  study)  for  the  integration  of  the  boundary 
layer  equations  and  the  process  for  obtaining  the  final  solutions 
for  the  shear  stress  and  velocity  profiles  at  each  location  along 
the  body  were  revised. 

In  the  sections  that  follow,  the  necessary  equations  and 
the  method  for  their  solution  are  developed.  The  relevant  computer 
program  presented  in  Appendix  A  was  executed  for  a  sphere  and  an 
ellipsoid,  and  the  results  are  discussed  in  Chapter  IV. 
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CHAPTER  II 

DEVELOPMENT  OF  EQUATIONS 

Boundary  Layer  Equations 

The  equations  to  be  solved  cover  the  case  of  axisymmetr i c. 
Steady  flow  past  a  blunt-nosed  body  of  revolution.  The  curvilinear 
coordinates  of  a  point  P  in  space  are  taken  as  (x,  y,  6).  The 
basic  notation  and  scheme  of  coordinates  is  shown  in  Figure  1;  U 

OO 

is  the  free  stream  velocity,  and  U (x)  is  the  velocity  in  the 
x-direction  just  outside  the  boundary  layer.  Theta  (6)  is  the 
angle  between  a  fixed  meridian  plane  and  the  meridian  plane 
containing  P.  The  surfaces  Tj  «  constant  and  T2  ■  constant  are 
taken  as  surfaces  of  revolution  about  the  axis  OX,  and  are  such 
that,  if  C  is  the  curve  of  intersection  of  the  surface  of  the  body 
by  a  meridan  plane,  then  the  sections  of  the  surfaces  Tj  »  constant 
and  f2  ■  constant  are  normals  to  C  and  parallel  curves  to  C, 
respectively.  Consequently,  T|  may  be  taken  as  the  distance  x  from 
the  forward  stagnation  point  0  measured  along  C,  and  as  the 
normal  distance  y  from  the  surface.  Here  u,  v,  and  w  are  the 
components  of  velocity  in  the  directions  of  increasing  x,  y,  and 
0,  respecti vel y ;  and  rQ,  the  body  radius,  is  the  distance  from 
to  CX.  so  that  r  is  a  function  of  x  alone.  The  coordinates 

'  ft 


(x,  y,  6)  form  a  mutually  orthogonal  system  with  the  following 


metric  coefficients: 

hx  “  1  +  Ky> 
h  =1 

y 

and 

he  "  r  ■  r0  *  V  cos  °> 

where  K  (a  function  of  x)  is  the  longitudinal  curvature  and  r  is  the 

radial  distance  from  P  toOx'*. 

For  flows  at  high  Reynolds  numbers,  around  bodies  whose 

local  radius  rQ  is  large  compared  to  the  boundary  layer  thickness 

6  (that  is,  6/rQ  <<:  !)  and  whose  surface  contains  no  large  variations 

2  2 

in  longitudinal  curvature  (e.g.,  sharp  corners  where  d  rQ/dx 

.  Q 

becomes  infinite)  so  that  dK/dx  -  1,  the  product  K6  is  small"'. 

Hence,  the  above  metric  coefficients  can  now  be  approximated  by: 

(1) 


boundary  layer  ’ 

=  0,  the  boundary  layer 

equations  are:  ’’ 


h  -  1  , 

y 

and 

he  3  ro> 

where  the  domain  of  y  is  restricted  to  the 

Since  the  motion  is  independent  of  6  and  w 
it.  10.  12 


u  ♦  v  *  U  4^-  +  V  4— *■  (Momentum)  > 
ox  dy  dx  >2 

7  ay 


l(ur  )  a (vr  ) 


0  (Conti nui ty) . 


The  boundary  conditions  are,  with  the  subscript  w  denoting  wall 


condi t i ons : 


y  *  0  :  u  *  0 
'  w 


v  *0 
w 

y  — ►  “  :  u  — ►  U  (x) 


Equations  (2)  and  (3)  can  be  combined  into  a  single  equation 
through  the  introduction  of  the  stream  function  Let 


u  ■  —  |—  (i pr  )  - 

rQ  3 y  T  o  dy 


)  3  \  _  .  M  i _ o 

7~  37  lvV  ^7  r  dx 

o  o 


The  resulting  equation  is: 


2£  it. .  (h  +  i_  tlo  )£±  .U  £  +  v  , 
ay  W  V3x  r0  dx/av2  dx  av3 


with  the  boundary  conditions: 


V  •  0  :  (fjw  ■  »• 


i  vm  m  •ncm  *•  y 


S' 


9 


f'"  -  -M.  (10) 

w 

Considering  Equat:on  (8),  when  R  -  0  and  H  is  constant,  the 

bracketed  expression  containing  the  partial  derivatives  with  respect 

to  x  vanishes  and  the  following  ordinary  non-linear  differential 

4 

equation  is  obtained  : 

f"*  „  .  +  M(f'2  -  1). 

This  relationship  provides  a  family  of  "similar"  solutions. 

A  final  transformat  ion  is  applied  to  Equation  (8)  for  the 
purpose  of  numerical  calculations.  The  quantity  4>  is  introduced, 
such  thct- 


and 


f  -  <P  +  n, 

f'  -<*>'  +  1,  01) 

f"  -  <t>",  etc. 


Equation  (8)  becomes: 

4>‘"  *  -  ♦  8^  (<>+n)  <p"  +  m (4> 1  ^  +  2ip' )  +  x  £(<p'  +1)  —  <j>"  -|^j, 

(12) 


with  the  boundary  conditions: 


and 


rt  ■  0  :  A  -  0,  <j>'  -  -1 

Tw  w 

n  — ►  00  :  $ 1  — *  0. 


(13) 


Equation  (12)  has  several  noteworthy  advantages  over  other 

possible  representati ons :  one  is  the  fact  that  the  starting  process 

is  very  simple.  At  x  ■  0,  all  x-dep?ndence  is  removed,  leaving 

only  an  ordinary  differential  equation  to  solve.  A  second  advantage 

of  Equation  (12)  is  that  almost  all  of  the  variation  in  boundary 

layer  thickness  has  beer,  eliminated.  The  thickness  in  the  transformed 

system  seldom  varies  by  more  than  50  percent  over  a  range  of  x, 

whereas  the  actual  physical  boundary  layer  thickness  might  vary  by 

a  factor  of  10  or  higher.  Finally,  it  is  very  important  to  note 

that  these  equations  are  entirely  free  of  mathematical  “pathologies." 

They  are  well  behaved  at  x  -  0,  and  solutions  of  Equation  (12) 

2  1 3 

exhibit  an  asymptotic  convergence  nature  for  large  g  ’ 


Displacement  Thickness  and  Skin  Friction 

1 2 

The  displacement  thickness  6j  is  given  by  : 

r 


1 


0  -  5)  dy. 


In  non dimensional  form,  this  becomes: 


<5 


* 

i 


(1 -f ' )  dn. 


o 

After  utilizing  Equations  (9)  and  (11),  the  nondimensi onal 
displacement  thickness  becomes: 


6  j  “  ! im  (n-f) 


The  local  skin  friction  coefficient  is  defined  as  the 
ratio  of  the  local  shear  stress  t  at  the  wall  to  the  local  dynamic 
pressure  outside  the  boundary  layer.  The  relation  is: 


2  )  f"  *  2  (■—)  6" 

llx  w  {Ux'  V 


where 


-  (y*  fw  -(4)4  *»• 


The  parameters  6^  and  f^  are  very  useful  in  studying  the 
growth  of  the  boundary  layer  and  the  variation  in  skin  friction 
along  a  body's  surface.  Graphs  and  some  additional  discussion  of 
these  quantities  for  a  spherical  and  an  ellipsoidal  body  are 
presented  in  Chapter  IV. 


Vorticity  Relationships 

The  coordinate  system  illustrated  in  Figure  1,  the  assump¬ 
tions  made  on  Pages  **  and  6,  and  the  metric  coefficients 
Equation  (1)  are  used  here  to  develop  the  vorticity  equation.  It 
is  assumed  that  v  is  negligibly  small,  and  the  magnitude  of  the 
gradients  in  the  y-direction  is  much  greater  than  the  x-wise 
gradients.  Also,  since  w  -  0  and  the  0-wise  gradients  are  zero 
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If  this  term  is  rewritten  as: 


u-u  o  (2nr  ) 
b  o 

2TTr  dx  '  * 
o 

then  it  clearly  represents  a  change  in  vorticity  due  to  the  extension 

or  contraction  of  vortex-lines  resulting  from  the  variation  in  body 

circumference  with  x.  There  is  an  intensification  of  vorticity  due 

to  the  extension  of  vortex- lines  when  d(2rrr  )/dx  is  positive*. 

o 

Consequently,  the  vortex  stretching  term  is  a  source  of  vorticity 
for  bodies  whose  radius  increases  downstream. 

In  order  to  study  the  terms  of  Equation  (19),  it  is  again 
convenient  to  introduce  the  stream  function  ^  for  u  and  v  and  the 
transformation  defined  by  Equation  (7).  The  non d i mens iona 1 i zed 
terms  of  the  vorticity  transport  equation  become 


i 


J 


(1+K+2R) 


0 


f"‘ 


(u-convection) , 


(v-convect i on) , 


V  dr0 


(stretch i ng) , 
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CHAPTER  I  I  I 
METHOD  OF  SOLUTION 


Basic  Scheme 

The  method  of  solution  that  is  used  in  this  work  is  based 
on  a  method  developed  by  Smith  and  Clutter'^.  They  replaced  the 
x-wise  partial  derivatives  in  Equation  (12)  by  finite  difference 
approximations  as  originally  suggested  by  Hartree  and  Uomersley 
(1937)-  The  equation  is  converted  to  an  ordinary  differential 
equation.  If  the  flow  in  the  x*direction  is  divided  into  n 
stations,  then  the  ordinary  differential  equation  has  to  be  solved 
n  times  in  succession,  since  M  and  R  as  well  as  <J> ,  and  $>"  vary 
with  x.  The  numerical  methods  used  are  a  form  of  the  ordinary 
finite  difference  treatment  since  discrete  variable  approximations 
are  made  in  both  the  x  and  n  directions. 

Solution  in  the  x-Direction 

Two  finite  difference  representations  of  Equation  (I?)  are 
possible.  The  first  is  called  the  "point"  form.  Here,  the 
differential  equation  is  written  to  apply  at  a  point;  that  is,  the 
x-der i vat i ves  at  a  point  are  replaced  by  their  finite  difference 
equivalents.  The  second  treatment  is  to  deal  in  terms  of  mean 
values  of  the  variable  0  over  a  finite  region;  this  is  referred  to 
as  the  "mean"  form. 
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After  extensi1  :  calculations,  using  both  the  "point"  and 

"mean"  forms,  Smith  and  Clutter  found  that  tne  "point"  forms  were 

generally  more  accurate  and  exhibited  a  greater  degree  of  stability. 

In  particular,  the  "point"  form  with  three  ana  four  points  proved 

1 3 

to  be  the  most  accurate  of  all 

The  basic  scheme  of  the  finite  difference  representat ion  is 
illustrated  in  Figure  2.  Tne  variables  M  and  R  are  known  as 
functions  of  x.  The  lines  ,  *n_].  *n-2’  etc>»  partition  the 
x,  n  space  into  a  number  of  regions.  Since  Equation  (!2)  is 
parabolic  in  form,  it  must  be  solved  in  the  direction  of  increasing 
x.  It  is  assumed  that  the  solution  $  (n)  and  its  derivatives  <p ' 
and  <J>"  (n)  are  known  at  all  stations  up  to  and  ir  '  iding  xn_}. 

The  problem  is  to  obtain  the  solution  at  x^.  To  accomplish  this, 
the  "point"  method  is  applied  to  Equation  (12)  at  x^. 

The  two,  three,  and  four-point  "point"  forms  that  replace 
thi  x-wise  partial  derivatives  in  Equation  (12)  are  easily  obtained 
by  differentiating  the  Lagranaian  interpolating  polynomials  of  two, 
three,  and  four  points,  respectively.  'n  the  actual  computer 
program,  these  formulas  are  written  tc  handle  unequal  step  lengths. 
ror  the  sake  of  simplicity  only,  the  four-point  formula  for  equal 
step  length  is  presented  here. 

The  four-point  formulas  for  n  and  ^)n  including  the 


error  terms  art: 
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x  )  n 


11$  *184  .  ♦  S$  ,-241  , 

n _  n-  1 _ n-2  n-3 

6Ax 


(Ax) \  U) 


(21) 


3x 


and 


lil)  «  ]l*il  l8^ri-l  *  9?n-2  '  2<J,n-j  +  (Ax)4  aV  (C) 


3x  /  n 


6Ax 


(22) 


3x 


where  x  ,  <  £(x)  <  x  The  error  terms  indicate  that  these 
n- j  —  —  n 

relationships  are  of  third-order  accuracy  and  become  exact  for 
third-degree  polynomial  variations. 

If  Equations  (21)  and  (22),  without  the  error  terms,  are 
substituted  into  Equation  (12),  the  following  equation  for  the 
four-point  approximation  results: 

C  ■  v^r  ♦  V  4  "n  K2  4  2*;>  *  ‘UK 

-  4  K,2  -  2K-l)  -  K  <"♦,  -  '«Vl  4  9V2  ' 

(23) 


Consequently,  Equation  (12)  has  been  transformed  into  an  ordinary 
non-linear  differential  equation  at  x^. 

It  is  important  to  observe  that  the  ratic  x/Ax,  rather  than 
the  step  length  Ax,  is  a  primary  parameter.  In  fact,  Smith  and 
Clutter  found  that  this  ratio  should  not  exceed  25  for  a  stable 
solution  using  single-precision  arithmetic  on  a  digital  computer. 


r  f '  rP#F5*W,?!;  *>,»  yKWMW* 


I**: 
Jt 

M 

$ 

£ 


)9 


Solution  of  the  Ordinary  Differentia)  Equation 

The  ordinary  differential  Equation  (23)  is  solved  with  a 
predictor-corrector  (or  ext rapoi at i on- i nte rpol at i on)  technique. 

2 

The  formula*  used  in  this  method  are  developed  in  Collatz1  book  . 
The  four-point  predictor  relationships  written  in  Lagrangian 
form,  including  the  local  error  terms,  are: 

♦«..  -  ♦,  *  h*i  *  r  *'i  *  75J  [,S8*V  -  'H»m  *  n*'i-2  -  Wi-j] 

*  0(h7), 

*|  +  ,  “  <?;  +  +  3So  [32^'iM  *  26^V-l  *  ,5S*iV  38*i'-3l  +  0(h6)* 


(25) 


and 


$'i‘+l  -  0\'  *  jz[h<p\“  -  59^:,  +  37*v:2  -  9*y:3]+  o^5)-  (26) 

The  correspond i ng  corrector  type  equations  are: 

2  j  p 

°i  +  1  "  *i  +  h0i  +  J~  0'{  *  L’7^1  +  ,20$'i”  ‘  21  ^i- 1  +  Hyi^th7), 


(27) 


♦i.i  -  ♦;  *  h*y  *  5S0  [36»'i';i  ♦  I7i*r  -  36»v-i  ♦  wv:J- oih6)’ 


and 


♦V.  1  -  *'1'  *  >9 -  5 *  *y:J  •  o«.*>. 


(28) 


(29) 
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In  tnese  formulas,  h  is  the  step  length  in  the  n-direction  and  i 
is  the  step  count.  The  procedure  is  to  make  a  prediction  by  means 
of  Equations  (24),  (2£)  ,  and  (?.6).  Using  these  values  in  Equation 
(23),  is  computed;  then,  Equations  (27),  (28),  and  (29)  are 

employed  to  obtain  improved  results.  As  a  consequence  of  the  fact 
that  the  above  relationships  provide  a  high  degree  of  accuracy, 
this  single  improvement  is  sufficient  for  each  n  step. 

The  solution  is  started  at  the  first  point  by  a  Maclaurin 
expansion  which  is  used  as  a  predictor.  The  ser<es  for  p  has  the 
form; 


<f>  -  A 


ny 


h2 

P" 

2  w 


h3 

rC 


(30 


and  similarly  for  p1  and  p".  The  values  of  4  ,  p'  and  p1’  are 
P  P  p  p  p 

substituted  into  Equation  (23)  to  find  the  extrapolated  quantity 
•Vp 1  •  These  extrapolated  values  are  used  to  compute  improved  values 
by  Obrechkoff's  corrector  type  formula  (Hi  1 deorand,  1956): 

2  2 

♦c  .  $w  +  I  %  *  V  ’  TO  (*p  -  V  +  120  %'  ♦C)  -  0(h?^  <31) 


with  corresponding  relations  for  p^  and  p^'.  An  imp'-oved  value  for 
p‘"  is  now  computed  from  Equation  (23).  The  improved  quantises 
and  p"  are  now  substituted  for  the  predicted  results  6  , 

WWW.  ■  p 

$p  and  lnt0  equations  of  the  type  (31)  to  produce  further 
Improved  values.  This  iterative  process  is  repeated  several  times 
until  the  convergence  of  p1^’  is  effected.  This  method  of  integration 
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Is  now  extended  to  the  next  two  points  by  using  two-  and  three- 
point  formulas  of  the  same  family  as  the  four-point  relations 
presented  above.  Finally,  the  four-point  formulas  can  be  utilized. 

A  constant  r,-step  iength  h  is  used  except  for  the  first  four 
points  where  its  values  is  b/k.  The  same  n-step  sizes  are  used 
throughout  the  entire  x-range  since  most  of  the  variation  in 
boundary  layer  thickness  has  been  removed  (cf.  Page  10). 

Integration  Procedure 

In  general ,  a  non-linear  boundary  value  problem  with  two- 
point  boundary  conditions,  one  point  being  an  asymptotic  condition 
at  n  ■  <»,  can  be  rather  difficult  to  solve.  However,  in  this  case, 
a  simple  method  consists  ->f  solving  it  as  an  initial  value  problem 
using  an  arbitrary  value  of  in?  unknown  condition  at  the  initial 
po.  it.  It  is  then  necessary  to  search  among  the  single  infinity  of 
possible  values  of  the  unknown  condition  for  the  particular  value 
which  allows  the  boundary  conditions  to  be  satisfied  asymptoti cal ly 
at  the  second  point. 

The  essence  of  the  present  method  is  to  seek  in  a  systematic 
manner  the  value  of  such  that  Equation  (23)  is  satisfied  in 
conjunct' on  wi th 

$■0,  -  - 1  at  T)  ■  0, 

and 


|VJ  <  e  at  n  «  n^, 
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where  e  is  a  prescribed  small  quantity.  The  detailed  process, 
illustrated  in  Figure  3»  is  as  follows; 

1)  Choose  an  initial  estimate  for  1“  and  a  value 

w 

for  n^. 

2)  Solve  Equation  (23)  by  "marching"  outward  away 
from  the  wall  using  the  technique  outlined  in  the 
preceding  section. 

3)  Cont i nue  unt 1 1 

(a)  >  -e,  or 

(b)  n  -  V 

k)  If  3(b)  occurs  before  3(a),  then  the  initial  estimate 
for  ^  was  too  low.  Therefore,  ^  is  increased  by 
an  appropriate  amount  a (  and  the  process  is  started 
over  again  at  step  2.  This  case  is  indicated  by 
curve  (a)  in  Figure  3> 

5)  If  3(a)  occurs,  then  there  are  three  possibilities: 
(a)  /$' /  <  e  and  n  •  When  this  occurs,  an 
acceptable  solution  for  ^  has  been  found. 

This  situation  is  represented  by  curve  (b)  in 
Figure  3. 


(b) 


>  e  (high  trial).  This  indicates 
initial  estimate  of  was  high.  It 
by  a  small  amount  a^  and  the  process 
again  from  step  2.  See  curve  (c)  in 


that  the 
is  reduced 
is  repeated 
Figure  3- 


2^ 


(c)  <  e.  Continue  "marching"  until  either  5(a) 

or  5(b)  occurs  or  becomes  less  than  -e.  If 
<  -e  (low  trial),  then  the  initial  value  for 
was  low  and  a  new  estimate  is  found  by  adding 
a  small  increment  and  then  the  procedure  is 
started  over  again  at  step  2  (curve  (d)  in  Figure 

3). 

6)  When  a  high  and  a  low  trial  are  found,  a  new  value 

for  d>"  is  obtained  from  the  arithmetic  mean  of 
w 

these  trials.  This  averaging  will  lead  to  conver¬ 
gence  because  varies  almost  linearly  with  <t>" 
for  a  wide  variety  of  flow  problems  investigated  by 
Smith  and  Clutter. 

7)  After  step  6  has  been  executed,  the  estimates  for 

are  continually  improved  by  taking  the  mean 
of  the  best  high  and  low  pairs. 

8)  When  three  values  of  $"  are  obtained  that  satisfy 

w 

the  condition  that  I  <p'  I  <  e  and  n  “  nw  then  the 
final  solution  is  found  by  forming  a  simple  three- 
point  interpolation. 

9)  Proceed  to  the  next  x-station  and  use  the  final  value 
of  from  step  8  of  the  previous  station  as  the 
initial  guess  for 
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It  should  be  mentioned  that  this  method  Is  not  completely 
devoid  of  difficulties.  As  the  Integration  proceeds  downstream 
in  the  x-dlrection,  a  great  sensitivity  to  initial  estimates  of 
<J>^  develops,  in  fact,  it  becomes  virtually  impossible  to  find  any 
solution  where  is  kept  within  the  range  +  e.  This  problem  is 
most  probably  a  result  of  the  following  causes:  0)  the  value  of 
and  the  associated  asymptotic  boundary  condition;  (2)  the  value 
of  e;  (3)  the  accumulative  error  due  to  the  product  of  the  four- 
point-  i nterpolati on  errors  with  the  ratio  x/Ax;  (A)  round-off  error 
(5)  a  coupling  effect  of  all  or  some  of  these  factors. 

After  a  considerable  number  of  trial  calculations,  it  was 

found  that  this  difficulty  can  be  effectively  overcome  by  lowering 

and  increasing  e  by  small  increments-  However,  these  changes  in 

either  n  or  e  should  be  kept  as  small  as  possible  in  order  to 

minimize  the  attendant  error  in  the  final  value  of  A". 

w 


CHAPTER  IV 


RESULTS  AND  CONCLUSIONS 


Discussion  of  Results 

A  computer  program  was  written  to  solve  the  equation  of 
motion  (12)  with  the  associated  boundary  conditions  (13)-  The 
program  was  executed  for  the  cases  of  a  sphere  and  an  ellipsoid 
with  an  eccentricity  of  1/2  (semi-major  axis  a  -  2  and  semi-minor 
axis  b  »  1.73205)  in  a  meridian  plane.  The  nondimensiona)  displace¬ 
ment  thickness  6*  and  rj  [see  Equation  (16)J  for  both  bodies  of 
revolution  are  illustrated  in  Figures  b  through  7.  In  Table  I,  a 
comparison  between  the  present  method  and  that  due  to  Smith  and 
Clutter  is  given  for  the  sphere. 

In  these  graphs,  the  curves  denoted  by  R  ■  0  represent  the 
spectal  case  when  rQ  is  constant.  The  increase  in  boundary  layer 
displacement  thickness  is  a  consequence  of  the  fact  that  the 
vorticity  generated  on  the  body's  surface  is  transported  chiefly 
due  to  the  action  of  viscous  diffusion,  while  the  effect  of  convec¬ 
tion  towards  the  wall  Is  negligible.  Whereas,  in  flows  where  r 

o 

varies,  there  are  the  opposing  actions  of  convection  toward  the  wall 

and  viscous  diffusion  away  from  it.  Consequently,  the  boundary 

layer  grows  more  rapidly  and  the  skin  friction  (shear  stress  at  the 

wall)  is  smaller  in  magnitude  for  the  case  when  r  is  constant 

o 

compared  to  the  case  when  rQ  varies. 
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TABLE  I 


COMPARISON  OF  VALUES  OF  f"  ON  A  SPHERE  AS 

w 


CALCULATED  BY  THE  PRESENT  METHOD  AND  THE 
METHOD  DUE  TO  SMITH  AND  CLUTTER 


6° 

f" 

w 

Present  Method 

f" 

w 

Smi th  and  C  lutter 

0 

1 . 3 1 3^8 

VO 

OO 

VO 

30 

1  .26038 

1 .25888 

60 

1.07652 

1 .09011 

90 

0.63127 

0.6562 

100 

0.3^365 

0.3580 

104a 

105. 93 

0 

0 

Extrapolated  values. 
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In  Figures  8  through  1 5  •  the  o' i  mens  ion  less  terms  (20)  of 
the  vorticity  equation  (ip)  have  been  plotted  for  several  longitu¬ 
dinal  locations  on  the  sphere  and  ellipsoid.  Here,  the  effects  of 
each  term  can  be  observed  and  the  most  important  of  these  is 
described  below. 

In  the  region  near  the  forward  stagnation  point,  the  vorticity 
produced  due  to  stretching  is  almost  exactly  equal  in  magnitude  to 
the  u-convection  of  vorticity  as  indicated  in  Figures  8  and  12. 

This  is  a  result  of  the  fact  that  near  the  stagnation  point  u  and 
u)q  can  be  written  as: 


u  -  cxf 1 (n) , 


3u  ~  c 

“e  *  "  7?  "  ’  I xf  * 


where  c  is  a  constant  .  In  addition,  since  r  -  x. 


dro  -  ^6 


r  dx  x  ’ 
o 


and,  therefore, 


UW.  dr„  ,2v 

_ 0  O  ~  C  X  r  I  f  II 

r~  a; - r f  f  - 

o 


fc>B  -e2* 

- r f  f  • 
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Figure  8  Dimensionless  Terms  of  the  Vorticity  Equation  for  a  Sphere 
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menslonless  Terms  of  the  Vortlclty  Equation  for  a  Sphere 
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- U-CONVECTION 

-  V- CONVECTION 


x=  0.1734550 
e*  0  5 


Figure  12  Dimensionless  Terms  of  the  Vorticity  Equation  for  an  Ellipsoid 
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Figure  13  Dimensionless  Terms  of  the  Vorticity  Equation  for  an  Ellipse 


U- CONVECTION 
V-CONVECTION 
DIFFUSION 
STRETCHING  / 
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Here,  the  diffusion  of  the  vorticity  produced  in  lower  boundary 
layer  is  effectively  opposed  by  the  action  of  -j.v.  ction  normal  to 
the  surface. 

Further  downstream,  the  counteracting  effect  of  diffusion 
becomes  very  evident  (Figures  10  and  lM-  Since 


]_  /St  1  =  1 

p  I  "dy  J  w  =  p  "Sx  * 


from  Equations  (2)  and  (15),  it  is  easily  seen  that 


■  (P) 

\  dy  /  w 


So,  from  Equation  (32),  the  integrated  diffusion  term  is  proport:onal 
to  the  pressure  gradient.  Now,  since  3p/3x  is  negative  due  to  the 
acceleration  of  the  external  stream,  it  is  apparent  that  diffusion 
is  becoming  predominant  and  the  vortic'ty  production  in  the  lower 
boundary  layer  (stretching)  is  transported  away  entirely  by  diffu¬ 
sion  to  higher  levels. 
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Cortclus  Ions 

The  principal  results  and  conclusions  obtained  from  this 
study  are  summarized  by  the  following  statements: 

1)  A  computer  program  was  successfully  developed  to 
solve  the  incompressible  laminar  boundary  layer 
equations  for  the  case  of  steady  flow  at  high 
Reynolds  number  about  blunt-nosed  bodies  of 
revolution  where  there  are  no  large  variations  in 
longitudinal  curvature.  The  method  appears  to 
provide  sufficient  numerics!  accuracy  for  most 
engineering  applications  except  near  flow  separation 
where  the  accuracy  fails  to  about  two  decimal  places 
(cf.  Table  1).  Lastly,  the  program  is  capable  of 
providing  solutions  rapidly,  since  the  total 
execution  time  was  under  three  minutes  for  either 
the  sphere  or  the  ellipsoid. 

2)  The  local  tendency  of  large  vorticity  production  in 
the  lower  boundary  layer  creates  larger  curvature  of 
the  vorticity  distribution  in  the  middle  boundary 
layer.  Therefore,  there  is  more  rapid  diffusion 
which  increases  the  boundary  yer  growth  and 
counteracts  the  trend  towards  large  values  of 
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vorticity  at  the  wall.  Hence,  Og  is  increased 
by  only  a  small  percentage,  even  though  the  local 
stretching  effect  can  be  very  large. 

3)  Vortex  stretching  tends  to  produce  a  greater  increase 
in  vorticity  in  a  region  close  to  the  surface. 
Consequently,  skin  friction  is  increased,  creating 
more  rapid  growth  of  the  boundary  layer  which,  in  turn, 
tends  to  inhibit  the  increase  in  skin  friction. 
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COMPUTER  PROGRAM  FOR  SOLVING  THE  BOUNDARY  LAYER  EQUATIONS 

The  following  computer  program  was  written  in  the  Fortran  IV 
G  language  for  the  IBM  360/67  digital  computer  of  the  Computation 
Center  at  The  Pennsylvania  State  University.  A  detailed  discussion 
of  the  method  used  can  be  found  in  Chapter  III. 


*5 


C  COMPUTER  PROGRAM  FOR  SOLVING  THE  BOUNDARY  LAYER  EOUATIUNS 

real  k 

01  MENS  ION  XI { 505) , X2 ( 5051 ,X3( 505 )  ,X4| 505 ) , X5( 3)  , SOLI! 3 ) , SOL 2 t  3 ) , $U 
1 L  3 ( 3 ) ,PREVO( 25, 505) , PREY  1 ( 25.505), LWl  300  1 , PI  I  1 00 ) ,Pp I ( 100 ) , PO I ( 100 
2 )  ,  P  T I  (  1 00 ) «  X  Y Z  (  300 ),EPS1(25), EPS 2(25), ST( 400) ,UV(25),TA(3),TB(3),V 
3SL(505),FSLN( 505) 

INTEGER  HH ( 450 ) 

DATA  X1/5#0./,X2/5*0./.X3/5*0./,X4/5*0./ 

PIA( A,P,C.O,E,F,G,H)  =  4*18-1. !-C*(0-E)+F*(G*H) 

PP  IA ( A.R.C.O, E, E )  =  -1. +A*( B*C ) -D* (E-F) 

PDI A ( A ,R ,C,0  )  =  A*R*(C*0) 

RFTA1(A,R.C.D,E,F,G,X,Y,2)  «  A*B*C +D* E+ F* { 1 88 . *G- 12 3 . * X+ 72 . * Y- 1 7. * 
1Z  I 

BETA2( A,8,C,0.E,F,G,X)  =  A*B*C* D* ( 3? 3 . *E- 264. *F* 1 59.  *G-38 . *X ) 
HFTA3(A,R.C  ,0,E,F)  =■  A+B*( 55.*C-59.*n+37.*E-9.*F > 

ZETA1  (A,R,C«D«E,F,G,X,Y,2  )  =  A+8*C+D*E+F* (  1  7.  *G*1  20  .  *X-2 1 .  *  Y-t-4.  *Z  ) 
ZETA21 A,B,C,0,6.F,G,X)  -  A+B*C*D* < 3R . *E+ 171 .* F-36. *G+ 7 ,*X ) 
ZETA3(A.B,C,0,E,F)  *  A*B*(9.*C  +  19.*f)-5.*E+F) 

RE AO  500.K.PHIDPW, ETAlNf ,H,P, Al, S 1 , L I M , S TUP  1 0 
RE An  501, ( EPS1 (HZ ) ,  HZ  *  l » L  I M ) 

REAO  501 , ( FPS21NZ I ,NZ=1,LIM) 

RE An  502, [CK.CK1,CK2,MESS,L0GM,C0GN,L0G0 
REAO  505,  (IIV1JL  ),  JL»1,LIM) 

REAO  506,1 rA(LAT),LAT*l,3),<T8(LBT),LBT*l,3) 

*2(  •  )  =  -1  . 

X  3  1  1  )  =  PMIOPk, 

FINANS  =  x  3 ( 1  ) 

XYZ< l )  »  10000. 

FSLN(l)  =  0. 

LSO  »  0 
LC  =  0 
HI  »  H/4. 

E  T  A  1  =  HI 
H 1  HAL  F  =  . 5  *  H 1 
H1S0  '  H ; H AL  F*M 1 
H1S01  =  . 2  *H ISO 
HI  CUB  =  HI *H1S01 / 1 2. 

HO  —  H/24. 

HSO  *  ,5*H*H 
HCUR  =  HSO*H/ 360 . 

HSO 1  =  HSC/180. 
no  24  LAM  x  l.LlM 
ICOUNT  =  0 

CALL  IDIOT (X, P, TA, T8,LAM,LOGM,LOGN,LOGO) 

I F ( LAM-MESS)  23.900,900 
900  LC  =  1 
23  LS  »  0 
IH  =  1 
IL  =  1 
MIN  *  1 
IX  *  l 
MST  =  1 

X4(  1 )  =  -EPS  1 (LAM) 

STU)  =•  X3(  1  I 

HH( 1 )  *  0 

LMU)  *  0 

11  PI  =  -HI+H1S0*X3( 1 )+(Hl»HlS0/3.)*X4U ) 

PPR 1  *  -l.4Hl*X3l  1 |4HlSO*X4! I  I 
P0PR1  =  X3U  !*H1*X4(  1  ) 

PTPRl  =  TR I  PI PREVO.FREV 1.LAM.EPS l (LAMI.ETAl ,2, PI, PPR 1 , P0PR1 ,EPS2(L 
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1AM j  .LC.tiv  I 

P[(l)  =  P 1A1H1HALF, PPR1 ,H1SCI,PDPR 1, X  3 (  1)  .H1CU8, PTPK1 , X4( 1 )) 

PP1  (1)  =  PPI A(H1HALF,PDPR1 ,X3( 1 ) .HlSOl.PTPRl,  X4<  1  )  ) 

POJ(t)  *  P0IA{X3(l>,HlHALF,PTPRl,X4(t >) 

PT1  (1  )  *  TR  I  PIPREVO,  PREVWLAM.EPS!  (LAM)  ,  ETA  1,2,  PI  (  1  1  .PPI  (  1  )  ,POI  (  1  ) 
1,FPS2(LAM),LC,UV) 

1  COUMT  =  1  COUMT  +■  l 
I  «  2 

1  J  *  1-1 

PMl)  *  PIA(H1HALF.PPR1.HIS01.P0P«1,X3<  1  ),H1CUB,PTI(  J)  ,XM1  n 
PPM  I  >  =  PPU(HlHALF,P|.‘PRl,X3i  1  1  .HlSU'l.PTI  (  J>tX4(l  >  ) 

POI(I)  =  PI)IA[X3(  1  )  .H1HALF.pt  I(  J),  X4(  1  I  ) 

PTI ( 1)  *  TRIP(PREVO, PftEVl, LAM, EPS 1 (LAM) ,ET A1 ,2,  PI  (  I  )  , PPI ( I)  ,Pl)l  (  I  ) 
1,EPS2(LAM),LC«UV) 

JF(A8S( PTI  t 1  )-PTI ( Jl ) .LE.CK1 )  GO  TO  2 
I  =  I  +  l 
GO  TO  1 

2  ETA  *  ETA1 

XI ( 2 )  *  Pl(I) 

X2 ( 2 )  *  PPI (  I  ) 

X3(  2  )  =  POI  (  I) 
x* ( 2 1  =  PTIIII 

CALL  STAMCH(PREV0,PREVUXlI2l,X2(2J,X3I2).X4(2),ETA,HltEPSl(LAM)  ,L 
IAM,XI(S),X2<S),X3(5),X4(-5),X4(<>)  ,  X4( 3 > ,  XI ( 4 ) , X 1 ( 3 > , X2< 4 ) , XZ I  3 ) , EPS 
22ILAM) »X4( 1 ),CK2,lC,UV,X3I3),X3<4) ,F$LN) 

NM  =  6 
FTA  =  H 
ITILT  *  0 

3  PTA  «  ETA+h 
KK  =  MM-1 
KX  »  NM-2 

K  Y  *  NM-3 
KZ  =  NM-4 

Y 1  -  BETA1 ( XI ( KK ) ,H, X2(KK ) .HSO.X3I  KK )  , HCUP. X4 { KK ) , X4 ( K X ) , X4 ( X Y I , X4 
1<KZ  )  ) 

Y2  *  BFTA2< X2 (KK ) ,H,X3(KK) ,HSOl, X4( KK ;,X4(KX) , X4( KY ) ,X4(KZ > ) 

Y3  3  BETA3 ( X3(KK ) ,H0, X4( KK I , X4(KX ) , X4( KY ) , X4IKZ I  I 

Y4  =  TRIPf  PR E VO, PREY 1.LAM.EPS 1( LAM ) .ETA.NM. Y1 »Y2»Y3»EpS2(LAM),lC,U 
IV) 

XI (NM)  *  ZETA1(X1(KK),M,X2(KK I .HSO, X3IKK ) , HCUB , 24 , X4 < K< I , X4 ( K X ) , X4 
1 (KY  )  ) 

X2(NM)  >  ZETA2(X2(KK),H,X3(KK),HS0l,Y4,X4(KK),X4(KX) ,X4(KY)  ) 

X3(MM)  i  ZETA3(  X3( KK) , HO, Y4, X4(<K )  , X4( KX ) , X4IKY )  ) 

X4(NM)  =  TR IP( PREVO.PREV1 , LAM, EPS  1 ( L A M ) , P T A , NM , X 1 ( NM ) , X2 ( NM ) , X 3 ( NM 
1 >,FPS2(LAM),LC,UV) 

1 F ( ITILT. FO. t )  GO  TO  32 
IF(X2(NM) .GT.-K)  GO  TO  31 
IF(ETA.LT.ETAINF)  GO  TO  930 

IX  =  IX+i 

X  YZ  (  IX  )  =  AB  S  (  X2  (  NM  )  ) 

1F(XYZ( IX) .GE.XYZI IX-1 ) I  GO  TO  6"20 
XL*  *  X3i  l  ) 

Mpg  *  MIN+l 

6020  X3( 1 )  »  X3< I >♦.025 
GO  TO  6 

930  IFJABSI X2INM) ).LE.l .  )  GO  TO  932 
X3(l)  »  .S*(HIGH-,X3(1)  ) 

GO  TO  777 
9  32  MM  m  HR *  1 
GO  TO  3 
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31  ITILT  =  1 

32  IF(*2INM>.G1.F)  GO  TO  4 
IF(X2<NM).LT.-K)  GO  TO  5 
!F(ETA.GF.FTA|NF)  GO  TO  8 
MM  -  MM*1 

GO  TO  3 

A  IH  i  lH*l 
HM( [HI  =  NM 

1FCHHI IHI.uF.hhi IH-1) )  GO  TO  41 
HIGH  s  X  3<  l  ) 

41  X  3 {  l  )  =  X  3  <  1 ) -S 1 
GO  TO  6 
S  IL  =  Il*l 
M[N  =  M I M  + 1 
L  W  (  I  L  )  =  M  M 

IF«LW<  IU.LF.LWt  IL-l  )  »  GO  TO  51 
XIW  =  X  3 ( 1  ) 

51  X  3 ( 1  I  =  X  3 1  1 ) *A  1 
5  I F ( X 3 (  l  i  )  301,66,66 
66  [F<  iH.GT. 1.  ANO.MIN.GT. 1  I  GO  TO  7 
GO  TO  1  1 

7  X  3 ( 1 )  =  .5«(hIGH+XLW) 

777  MST  =  M$T*1 

ST(MST)  *  X3(  1) 

I F  <  ST (MST ) .PO.STIMST-l  )  1  CALL  CHGET A  I E T A  I NF , X 3 1 1 ) , 523 , E25 , F 1  NANS , L 
ISO, STUDIO. ICK  ) 

GO  TO  1  l 

8  I.S  =  LS*l 

SOL  3( LS )  =  X 3 ( mm ) 

SOL 1( LSI  =  X 1 <  NM ) 

SO  L  2 ( L  S  )  =  X2(NM) 

X6(LS )  =  X3 ( 1 ) 

I F ( L  S . FO  ,  3  )  GO  TO  10 

IF ( S0L2ILS ) .GT.O.  )  GO  TO  9 

X3( 11  =  ,5»!X5(LS1+H1GH) 

GO  TO  1 1 

9  X3<  11  *  .5*( X5<  LSl+XLM  ) 

GO  TO  1  1 

10  I F (  SOL 2(1 ) — SOL 212))  9865,9864,9865 

9864  IF(SUL2U ) -S0L2I3)  1  9865,9866,9865 
9866  !F(  S0L2I 2  I-S0L2I 3 )  )  98  65,9868,9865 
9868  FINANS  -  X5<ll 

GO  TO  200 

9865  FINANS  =  OFC I j ( X 5, SOL 2 ) 

200  PRINT  2000 

PRINT  2100, ( SOL l (LLS).S0L2( LLS ) , SUL  31 LLS ) , X5 1 LLS ) ,LLS= 1 , 3  1 
PRINT  2?nO,cINANS»LAM,fcT41NF,L50 
X3( 1 )  =  FINANS 
PRFVOI LAM, 11=0. 

PRFVKLAM.l  )  =  -1. 

API  =  -Nl+HlS0*X3l lI+CHl*HlS0/3. »*X4(1( 

A PPR  1  =  -l.+Hl-»X3(l  I*HIS0*X4U  ) 

APOPR 1  r  X  3H l  ) +H 1  *  X4 (  1  ) 

APTPR1  =  TR IP( PREV0*PRE¥1 .LAM.FPS 1  I L AM ) , F T A  1 , 2 , A P 1 , A PPRl , APOPR 1 ,£P 

1  S2  i  L  AM  )  ,  l_C  ,  l»V  1 

PI ( 1  <  *  PIAIHIHALF,  APPR 1.H1S01, APOPR l ,X3C 1 > , Ml CUB , AP T PR 1 , X4 1  1  )  ) 
PPI(I)  =  PPI AIH1HALF, APDPRI ,X3( 1)  ,HIS01 , APTPR1 , X4(  1  ,  I 
pni(l)  =  PO  I  A!  X3<  1 1  ,HH»LF,  APTPRl ,  X4(  1  )  ) 

PTK1I  =  TRIP!  PREV0.PRFV1  .LAM.EPSULAM)  ,ETA1 ,2,P1  !  1  I  ,PP1  II  I  ,POI  «  1  ) 

l,FPS2(LAMi,cC,UVj 
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% 


L  «  2 

201  LJ  «  t-1 

PHD  «  PIA< HIHALF.APPR  UH1SQ1  , APDPR1  ,X3(  l > , H1CUB , PT I ( L J ) , X4 l 1 )) 
PP1IU  «  PPI AfHlHALF, APOPR1 , X3( 1 ) .H1S0I ,PTI (LJ ) ,X4( 1 ) J 
POIIL)  »  Pt)IMX3m,H5H4LF,PTHU),)(4(ljJ 

PT  1  (  L  )  *  TK I P<  PREVO, PREVl.LAM.EPSl I  L AN ) , ETA  1 , 2, PI ( d , PP 1( L I i PO l < L ) 
UFPS2UAMI  ,LC,UV) 

IF(A8S(PTI(L  l-PTJILJ) I.UE.CK1  !  GO  TO  202 
L  *  t  +  i 
GO  TO  201 

202  FTA  =  ETA1 

XI (2)  »  PI(L) 

FSLN12I  *  X  I ( 2  )  +ET  A 
X2<2)«  PP  I  CL  J 
X3(  2)  =  PO!  (U 
X4 (2)  «  °T I ( L  ) 

PRFV0ILAM.2)  *  P I  <  l  ) 

PREV 1 ( LAM, 2 1  =  PPI(L) 

CALL  STANCH!  PREVO.  PR6VUXU2)  ,X2(2),X3(?I.X4<2),ETA.Hl,fcPLl(LAM),L 
UM,Xl(5l,X2(5),X3<5],X4C5),X4(4),X4m.Xl(4),Xl(3).X2(4),X2<3),EPS 
22UAM),X4m,CK2,LC,UV,X3(3).X3(4!,ESLN) 

LN"  *  0 
FTA  »  H 

203  ETA  ■=  ETA+H 
LKK  «  LNM-1 
LKX  *  LNM-2 
LKY  *  LNM-3 
LKZ  *  LMM-4 

Yl  =  BETA1 (XIILKK) ,H, X2( LKK) . HSU.X3ILKK ) , HCIIB , X4 < L K<  >  , X4 ! LK X ) , X4 I L 
1KY),X4(LK2 ) ) 

Y2  =  BETA?(X2(LKK),H,X3(LKK ) , HSQ 1 , X4 ( L KK ),X4(LKX),X4(LKY),X4|LKZ)) 
Y 3  =  BFTA3I  X3ILKK  ),Ht),X4(LKK  )  ,X4(  LKX)  ,X4(LKY)  ;  X4(  L<Z  )  ) 

Y  4  x  TRIP (PREVO.PREVl.LAM.EPSH LAM)  , E T A , LNM , Y 1 , Y 2 , Y 3 , E PS2 ( L AM J  ,LC. 
HIV) 

XI ( LNM  )  =  ZETA1 (X 1 (LKK ) ,H, X2( LKK ) ,MSO, X3ILKK ) ,HCUB. Y4, X4( LKK ) , X4<  L 
1KX ) *  X4( LKY  )  ) 

FSLN(LNK)  =  XHLNMI+ETA 

X2(LNM)  =  ZFTA2(X2(LKK),HfX3(LKK), HS01 . Y4,X4(LKK),X4(LKX),X4(LKY)) 
PRFVO(LAM.lnM)  =  XKLNM) 

PREV1 ( LAM.LNM)  =  X2ILNM) 

X  3  {  L  )  *ZETA3(  X3(  LKK  )  ,  HO,  Y4  ,  X4  (  LKK  ),X4<LKX),X4(LKY1) 

X4(LHM)  x  TRIP( PREVO.PREVl.LAM, EPS  1 (L AM) ,ETA, LNM, xll LNM) ,X2 (LNM)  ,X 
13(LNM),FPS2(LAM),LC,UV) 

IF(FTA.GE.ETAINF)  GO  TO  204 
LNM  »  LNM+1 

GO  TO  203 

204  PRINT  504.X2(LNM)  ,  X3(LNM) 

PRIKT  4000.  ICOUNT,lJV(  LAM) 

OISPL  *  -Xl(LNM) 

00  400  IKW  x  l.LNM 
400  VEL (  I K  W )  -  X2(  I  KM ) ♦ 1 , 

PRINT  5000 

PRINT  5001 , ( \ -L ( IKZ ) , IKZ«1 ,LNM| 

PRINT  6002  »  0  I  SPL 
PRINT  5003 

PRINT  5004, ( X3( LKW) , LKM«1 , LNM ) 

PRINT  5005 

PRINT  5006, ( FSLNILZZ ) ,LZZ*1 ,LNM ) 

PRINT  4007 

PRINT  5008 « ( X4 ( L BG ) , L 8G *  1 , LNM ) 
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24  CONTINUE 
GO  TO  25 

301  PH INT  503, X3< 1 1,1AM 

25  STOP 

500  FORMAT!  7F10.7  ,  I3.F5.2  ) 

501  FORMAT! 5F1 2.8  ) 

50?  FORMAT!  HO,  ?F  12.3,  1 10, 31  3  i 

503  FORMAT! '0' . <A  N6GAT 1 V6  TRIAL  VALUE  FOR  WALL  SHEAR  =  '.F14.8,'  OCCO 
1  RRED  AT  law  =  • , 14 ) 

504  FORMAT!*  '.'AT  INFINITY  X2  =  *,F14.8,*  AND  X3  *  *,F14.R) 

505  FCV'MA T( 6F1 2.5 ) 

505  FORMAT(6F10.7 ) 

2000  FORMAT!  *0' ,7x, • SOL  t ' , 1 6X , • SOL  2  ' , 1 5X , ' SOL  3 • . I 2X , • X5  •  I 
2100  FORMAT!*  ' , 4! F1A.8, 5X  )  ) 

2200  FORMAT!  • 0  >F  1NANS  *  • , F 1 4 . 8 . 5 X , • FOR  LAM  *  * , 1 3 , 5X , • E T A  1  NF  =  *,F14 
1.7,SX, ' TOTAL  NO.  OF  CHANGES  IN  ETAINF  *  >,I4) 

4000  FORMAT! *0* ,' (COUNT  =  *,I4.5x,*FUK  TmETA  =  *,F12.fe,‘  DEGREES*) 

5000  FORMAT! *0* . *ThE  VELOCITY  PROFILE  IS*) 

5001  FORMAT!  •  •  .  1  1  P 1 1  .7  ) 

500?  FORMAT! *0* , «THE  NON  DIMENSIONAL  DISPLACEMENT  THICKNESS  IS',F12.7) 

5003  FORMAT! *0' . 'THE  VELOCITY  GRADIENT  PROFILE  IS') 

5004  FORMAT! •  *.IIF11.7J 

5005  FORMAT!  <0* . 'THE  SOLUTION  PROFILE  FOR  F  IS*) 

5006  FORMAT! •  '.IIF11.7) 

5007  FORMAT! • O' THE  PROFILE  FOR  F  TRIPLE  PR|M£  IS*) 

5008  FORMAT!  '  1  ,  11FU.7) 

FNO 

SUBROUTINE  STAMCHIF ,G, A A,  BB . CC . 00, E . HO , YM , I . R I O , R  II , R I  2 , R T II . X X 1 , Y 
1V1,A0,30.A1,B1,P,W,C,JJ,U,X33,X34,FS) 

DIMENSION  F ( 25, 505 ) , G (  25 , 505 ) , PCP I  50 ) , PCR ( 50 ) , PCD  I  50 ) , STC ! 50  I , 0 1 0 1 
150  I  ,01 1 ( 50  )  ,012! 50 ) ,0T 11 50) ,U( 25) ,FS(505) 

INTEGER  P 

SIGO(A,B.C,n,E,F,G,X!  =  A*B»OC*E  +  F*(  1 7. *G  +  1 20  .»X > 

SIG1 ( A.R.C.D, 6,F  )  =  A  +  3*C*D»( 33,«E+171 . *F  ) 

S I G  2  (  A ,  ft ,  C  ,  0 )  *  A+ft*(9.«C+19.*0l 

GAMO!  A.R.C.D,  E,  F,G,  X,  Y  !  *  A-Fft  *C+0*  E  ♦  F*  (  l  7  .  *G*  1  20.  *X- 2 1  .  •  Y  ) 

GAM> ( A.R.C.D, £,f,G>  =  A*B«C +U* ( 33 ♦ 1 7 1 . wF-36. *G ) 

GAM2I A.R.C.D, E)  *  A*B»1 9.»C*19.»0-5.*E ) 

OEL  TO ( A , ft. C,D,E . F,G, X, Y,2 )  *  a+B*C*D*E*f» ( 1 7.*G*t 20 .• X-2 1 . * Y  +  4 ,«Z ) 

OELTl (A,R.C.D,E.f,G,X J  =  A  +  B*C*D* ( 33 , *E* 1 7 1 . »F-36 .*G*7 . *X  ) 

0ELT2! A,B,C.D,E.F)  *  A+R*( 9  .  *C  ♦  1 9  .»D-5  .  *F-»F ) 

E  =  f*hO 
HOD  =  HO/24. 

HOSO  =  HD»HO 
H0S01  *  ,S*HOSO 
HS02  =*  HOSO  1/180. 

MOCUB  =  H0*HnS0/720. 

PhIP|  *  AA+HO»Rfl*HOS01»CC*HOCUB*! 183. *00-123. *W) 

PhPPI  s  RR+hO*CC*hS02* ( 323,»00-254. *w ) 

PHDPl  =  CC*HOO* I  55 .*0  )-59.*W) 

STPl  =  TRIPIF.G.l, YM , E , 3 , PH IPl, PHPPI , PHUP 1,R,JJ,U) 

PCP(l)  =  S I  GO! AA.HO.Bfl. HOS01 ,CC , MOCUB, STPl , OD ) 

PCR(l)  *  SIG1 ( BB,H0,CC,HSQ2 , STPl ,DD) 

PCO(l)  =  SIG2(CC. HOD, STPl, OD) 

STCI1)  =  TRIP! F.G. I .YM.E.3, PCP! 1  I ,PCRI  1 ) .PCOI l  ),K, JJ.O) 

5=2 

101  L  ■  R-l 

PCPIRl  =  S IGO! AA.HO.BB. HOSO! ,CC .HOCUB, STC! L ) ,00) 

PC«(RI  *  S IG1 I BB ,HO,LC ,HS«2, S>L( L ) ,0D) 

PCDK)  =  SIG2(CC.M!1D,STCCL),00) 


50 


STCiK)  =  TRIPrF,G,l,YH,E,3.PCP<K).PC«tK),PCOlK),R,  JJ,U) 
IF{A8S(STC(K)-STCIL)).CE.CI  Gil  Til  102 
K  *  K  ♦  1 
GO  TO  101 

102  FS(  3)  «  PCPOU+E 

e  *  e+ho 

00  *  PCP < K >+HO*PCR (<) +H0S01*PC0(K )+HOCUB*t 188.*STC ( K ) -123. *00+72 .* 
1W) 

01  =  PCR (K  l+HO*PCO(K )+HS02*<  323. *STC(K >-264. *00+189. *WI 
02  «  PCOf  K  J  ♦•tfrO*  (58.*STC(K)-59r*DD+37.*W) 

0T1  s  TR1P(F,G,1,YM,E,4,C0,01,02,H,JJ,I)> 

010(1)  *  GAMOI PCP(K ),HO,PCR( K) . HOSOl . PC D ( K ) , HOCUB , 0T1 , STC (K ) ,00) 
Oil(l)  *  GAM1(PCR(K),HO,PCD(KI,HS02.0T1,STC(K),OOI 
012(1)  *  GAH2( PCO(K ) .HOD, 0T1, STC(K)  .00) 

OTItl)  *  TRIP! F.G, 1 , YM,E,4,0I0< 1  I .01 1 ( 1 ) ,012 (1) .8. JJ,U> 

N  *  2 

103  P  *  N~1 

OIO(N)  *  GANG ( PCP(K ),HO,PCR<  K) .HOSOl. PCDIK > .HOCUR.OTI (  P I , STC ( K ) , DO 
1  ) 

0  1 1 ( N )  -  GAHl(PC«(<),HOtPC0(K).HS02,OTM P).STC(K) ,00) 

QI2 ( N )  =  GAH2I PCO(K) .HOO.OT I ( P) ,STC ( K ) ,0D) 

CTI(N)  *  MlP<F,G,I.YN.e,4,OtOIN),QniN),OI2(N),R,JJ,U> 

IF< AfiSIOT! (N)-OTI (P) J.LE.C)  GO  TO  104 
N  *  N+l 
GO  TO  103 

104  F S ( 4 )  *  0I0(N)+E 
F  *  E+HO 

«0  »  OIO‘N)+HO*OH(N)+MOSQ1»Q12(N)+Mncun*l 188.*OTi(NI-123.*STC(K>+ 
I72.«nr)-17.*H) 

« 1  *  or 1( N) +H0*012(N)+HS02*( 323.*0TI l N)-264.*STC( K)  + 1 89. *00-38 . »w ) 
R2  *  0(2 (N)+HOO«< 58.*OT I <NI-89.*STC IK ) +37. *00-9. *W ) 

RTl  =  TRlP(F,G,l,VM,E,5,ftO,Rl,R2.ft.JJ,U> 

K ( 0  *  06LT0(0I0(N) ,H0, 0 IltNI, HOSOl, 012(H), HOCUR.KT 1,0TI|N),STC(K)» 
100) 

ft  1 1  =  PEt-THOmm, HO, 012(H), HSQ2. RTl, OTHNI.STCIK). 00) 

812  =  OEL  T2(OI2(N),HOO,  8T1.0T  UN) ,  STC  IK  )  ,0D) 

RT11  *  TR  IP< F.G. I ,YM,E ,5,R IO.RI 1 ,RI2,R, JJ.U I 
F S ( 5 )  *  RIO+E 
XXI  =  OTI (N) 

YY1  *  STC(K) 

AO  *  OIO(N) 

BO  *  PCP(K) 

A 1  *  OIl(N) 

B1  =  PCR(K) 

X33  *  PCO(K) 

*34  «  0 1 2 ( N ) 

F(MI  «  RO 
F(  1,4)  «  AO 
F ( 1,5)  *  RIO 
G  f  1,3)  *  «l 
G  (  I  ,  4  )  *  A  1 
G (  1,8)  «  R I  1 
RETURN 
END 

SUBROUTINE  CHGETA( X, V , •«•, 2 , I , u, J ) 

*  1*1+1 

IF!  |  .r.T.  J  )  GO  TO  20 
X  *  X-ll 

*1  V  •  2 

ft 6  TURN  l 


51 


20  PRINT  ROOO 
ft F TURN  2 

8000  FORMAT ( 'O' , 'FAILURfc' ) 

F  NO 

SUBftOtjTI  NF  I  !J1  Ul  I  X.T  ,  A.H.C  ,M,  N.  JU> 

OIMFNSION  A  I  3  )  ,  ft  !  3  J 
IF(L-M)  1,2,2 

2  IFIL-NI  3,4,5 

3  X  *  All) 

Y  «  HU) 

GO  TO  1 

4  IF(L-JO)  5.6,6 

5  X  .  A(2) 

Y  =  8(2  ) 

GO  TO  1 

6  X  *  A< 3) 

Y  «  R( 3) 

1  RF  TURN 

FNO 

Function  oec IG( X, Y) 

OIMFNSION  X(3),Y(3),P(3I,S<3).IXP(3),IXS!3) 
I  J  f*  =  0 
IKS  =  0 
J  *  0 
HPOS  a  1. 

K  =  0 
BNFG  =  1. 

00  5  I  =  1,3 
IflYdll  3,5,1 

1  J  =  J*1 

P(  J)  *  Y(  I  1 

JXPIJJ  x  I 

IF( P( J) -HPOS )  2,5,5 

2  BP0S  =  P(J) 

:j  =  i x  p i j  i 
ijp  =  i 

GO  TO  5 

3  K  *  K*l 

S ( K )  X  ABS(Y(  I  )  ) 

I  XS(K I  *  I 
IF< S(K 1-ftNFG)  5,5,5 
5  BNFG  «  S(K) 

X K  =  ! XS ( K  ) 

IKS  =  1 

5  CON  T  I  N,  IF 

ifi  i jp.eo. n  go  to  5 

OFCIO  =  X(KK ) 

GO  TO  111 

6  I F ( lKS.FO.l ]  GO  TO  7 
OFC ID  *  X  <  J J  ) 

GO  TO  111 

7  IF  |  ABS< X( JJ)-X( KK ) J-. 0000001 )  1  3,  13,  77 
13  I F ( ft  POS-ftNFG I  131,131,133 

131  OFCIO  -  X( JJ 1 
GO  TO  111 
133  OFCIO  «  X(KK) 

GO  TO  1 1  l 

77  I F ( flPOS-BNEG  I  8,9,10 
HR*  BNFG/ftPOS 
RAT  «  A  INI ( R  I 


OECin  =  (RAT*X< JJ ) +  XIKK  '  1/ (RAT+1  .  ) 

GO  TO  111 

9  DFCfn  *  .  5*  ( X(JJ)*X(KK)  ) 

GO  TO  111 
10  It  =  BPOS/BNFG 
RAT  *  AINTIR  1 

OeCID  =  <RAT*X< KK ) +X< JJ) )/(RAT+l.) 

Ill  RETURN 

pno 

FUNCTION  TRIP(X,Y,J,XM.T,K,A.B.C.R,N.S) 

DIMENSION  X(25,  505  ), Y(25,505) ,S<  25) 

Z  =  -( . 5*( XM*1. )*R )»( AfT !*C+XM*I B*R+2.*B ) 

IFIJ.GT.3)  GO  TO  704 
GO  TO  ( 701. 702, 703), J 

701  TRIP  *  l 
GO  TO  705 

702  TRIP  -  2-M  (B*l.  )*(  B-Y<  1,K)  )-C»lA-X<  1  ,K  I  )  ) 

GO  TO  705 

703  FO  =■  $(  J)*(2.«Sf J) -S(  J-2)-S<  J-l)  )/(  (S(  J)-S(  J-2)  l*(S(JI-S(  J-1U) 

FI  =  SI Ji*lS( J  >  —  S ( J-2) )/i ( S( J-I )-S< J-2) )«( S< J-l )-S( J  )  )  ) 

F2  =  S( J 1* (S( J)-S< J-l ! ) / ( ( S< J-2)-S( J-l 1 )*! $ ( J-2 ) -S I J ) ) > 

TRIP  =  Z+ I B  +  l . I* (F0*B*F1*Y( 2,K) *F2«Y( 1 . K ) ) -C * ( FO* A+F 1 * X ( 2 , < ) +F 2 *X  i 
11. K)  ) 

GO  TO  705 

704  0  *  S  (  J  )*HS(  JI-SI  J-2  )  )*  (  S  (  J) -St  J-l  )  )) /(  <S<  J-3I-SI  J-2  I  )  *!  S  I  J-3  )-S< 
1 J- 1 ) )»( S  ( J-3 ) -S I J ! )  ) 

ll  =  S(J)*((S(J)-S!  J-3 ) >*IS( J  J -S I J-l ) ) ) /( (S(J-2)-S(J-3) I  *  (  S  I  J  -  2 )-S( 
1J-1) )*( S ( J— 2 ) -S ( J ) )  I 

V  «  S  (  J I  *  (  ( S ( J • -SI  J-3) )*IS(J)-S(J-2)))/((S(J-l)-S<J-3))*(S(j-l)-S( 
1  J-2 ) )»( S( J-l ) -S ( J ) )  ) 

H  *S« J)*(  (  S( JI-SI  J-3) )•( S( J )-S( J-2) )  ♦<  SI J )-S( J-3) >«<  S( J)-S( J-l  )  >  «•  1 
1 S ( JI-SI J-2 ) )*(S<  j )  — S f J-l ) ! i/( I  SI J i-S( J-3 I )*(S ( j  >-S<  J-2 )  l*( S( J  )-S( J 
2-1  ))  ) 

TRIP  =  Z*( 3*1 . )*<R*B*V*y ( J-l ,K i «U*Y( J-2,<) *0*Y( J-3.K ) ) -C» i w*A+ V*X ( 
l J-l ,K 1 +U*X ( J-2,K ) *Q*X ( J-3.K  )  I 

705  RETURN 
END 
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The  input  data  necessary  for  the  proper  execution  of  the 
above  program  should  consist  of  the  following: 

1)  The  first  data  card  should  contain  values  for  +£, 
an  initial  guess  for  f*J,  ri^,  the  n-step  size  (h) , 

-e,  a,,  a„.  the  number  of  x-steps  and  the  increment 
by  which  r;  is  lowered.  (See  pages  21-25) 

2)  The  next  set  of  data  cards  consists  of  M  and  R,  in 
that  order.  These  are  obtained  from  the  potential 
flow  and  body  geometry  for  the  particular  body  under 
consideration.  (See  Appendix  B) 

3)  The  next  data  card  consists  of  the  maximum  number 
of  allowable  changes  in  two  tolerances  for 
determining  numerical  equality  (1x10  ^  for  single 
precision),  an  approximation  for  the  number  of  the 
particular  x-station  which  is  nearest  to  the  point 
of  separation,  and  the  numbers  of  the  particular 
x-steps  where  the  values  for  +e  and  -e  are  to  be 
changed . 

4)  The  next  set  of  data  cards  cons’Sts  of  the  numerical 
values  of  x  (i.e.,  arc  length)  at  each  x-step. 

5)  The  final  data  card  should  give  the  numerical  values 
for  +c  and  -e  at  each  step  where  they  are  to  be 
changed  (cf.  statement  number  three  above). 
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APPENDIX  B 

CALCULATION  OP  M  AND  R  FOR  AN  ELLIPSOID 


Coordinate  System  and  Geometric  Relations 

A  semi -el  I ipt ic  coordinate  system  is  a  very  convenient 
representat i on  for  obtaining  the  potential  flow  solution  about  an 
ovary  ellipsoid.  The  defining  relations  are: 

X  *  Ky6, 

Y  «  K(l-u2)i  (52-l ) ^  cos  u,  (Bl) 

and 

2  -  K(l-L2)i  (62-))i  sin  u, 

where  -1  <  y  <  1  and  I  <  6  <  Tht  surfaces  5  ■  constant, 

y  =  constant  are  confocal  ellipsoids  and  hyperboloids  of  uo 
sheets,  respectively,  with  common  foci  at  (+K,  0,  0),  and  uj  is 
taken  to  be  the  azimuthal  angle  in  the  meridian  plane.  The 
coordinates  y,  6,  u  form  an  orthogonal  system  with  the  following 
metric  coefficients: 


and 


( B  2 } 


(B3) 


4 
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and 


h  ,  -  K(!-u2)i  {(]2-l)i  . 

CO 


(Bk) 


Considering  the  meridian  plane  with  Z«0  of  the  ellipsoid  <5 
with  semi-major  and  minor  axes  a,  b,  respectively,  the  following 
geometric  relationships  hold: 

.  _  (a2-b2)* 


(B5) 


and 


K  -  ae, 


(B6) 


where  e  Is  the  eccentricity.  The  equation  of  the  ellipse  in  the 
meridian  plane  is  given  by  the  familiar  equation: 


X2  Y2 

2.2“ 
a  b 


(B7) 


Now,  from  the  relations  (Bl)  with  (85)  ,  (B6) ,  and  (B7) ,  it  is 
easily  proven  that  for  this  body  (Lamb,  19^5): 


and 


6  ■  — 

o  e 


X 

y  *  7  • 


<B8) 


(B9) 
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Calculation  of  R 

Referring  to  Chapter  II,  the  expression  for  R  can  be 
..'it ten  as: 

p  x  dro  dX 

Ra7~dT  d7  • 

o 


(810) 


Now,  from  we  1 1 -known  formulas  of  Integral  calculus,  the  differential 
of  arc  dx  and  the  arc  length  x  are  given  by: 

i 


dX 


(fill) 


and 


x  »  a  E  (e,y) , 


(B12) 


where  E  (e,y)  is  an  elliptic  integral  of  the  second  kind  defined 
by: 

i 


E  {( 


-J'fw 


dy. 


(813) 


dr 


JV 

Since  rQ  «  Y,  is  readily  obtained  from  (B7)  and  gj-  from  ( B 1  1 ) 
Finally,  the  equation  for  R  becomes: 


-  y£  (e,y)  £{l-y2) (l-e2y2)J 


(814) 
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Calculation  of  H 


The  velocity  pctemia’i  4>  for  an  ellipsoid  moving  parallel 
to  the  X  axis  is  given  (see  Lamb,  19^5)  by: 

<?  «  Ay  £i6  in  jq-  -  1J  l 


where  A  *  all 


[i? -i 


From  (BIS)  and  (B2)  ,  the  velocity  in  the  y  direction  V  is: 


,r  .  5+1  , 

In  TT  ’  1 


For  the  ellipsoid  6  ■*  constant  -  6q,  the  relation  for  « 


takes  on  the  following  form: 


The  quantity  -s^-  on  is  obtained  from  (B16)  and  is  found  to  be: 

G/)6  ■  ?  °-'2)  U  [('VlO-cV)3]  1  [i  jl|  -  ) J.  <Blf 

o 

The  expression  for  on  is  most  conveniently  established  by 
differentiating  the  coordinate  relations  ( B 1 )  and  solving  the 
resulting  set  of  simultaneous  equations  algebraically.  The  result 


J 


(BIS) 


After  making  the  appropriate  substitutions  into  (B  5  7) .  M  finally 
becomes : 


M=  (1-e2)  yE  (e ,y)£(i-y2)  ( 1  -eV) 3 ] 


(B20) 
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